void resid_one(double *res, double *u, double *rhs, int n)
//Returns minus the residual for the model problem. Input quantities are u[1..n] and rhs[1..n], while res[1..n]is returned.
{
  int i;
  double h,h2i;

  h=1.0/(n-1);
  h2i=1/(h*h);
  //interior points
  for (i=2;i<n;i++)
    res[i]=-h2i*(u[i-1]+u[i+1]-2*u[i])+rhs[i];
  //boundary points
  res[1]=-h2i*(u[2]-2*u[1])+rhs[1];
  res[n]=-h2i*(u[n-1]-2*u[n])+rhs[n];
}
